!-----------------------------------------
!  A. Gazizov, LNGS, 18.12.2010
!  included in project ZN on 26.11.12, LNGS
!  $\gamma + p \to e^+ e^- p is added
!-----------------------------------------
module approx1x
   implicit none
   integer, parameter:: ne  = 264  &
                      , nm  = 3    &   ! <(1-x)^m>, m = 0, 1, 2
                      , nc  = 4    &   ! pp, pn, nn, np
                      , nee = 314

   real(8), dimension(ne)       :: egg
   real(8), dimension(nm,nc,ne) :: P
   real(8), dimension(2)        :: lnegg, lneep, lnavreep
   real(8), dimension(nm,nc,2)  :: lnP
   real(8), dimension(nee)      :: eep, avr1xeep

!------------------------
         CONTAINS
!------------------------

!-----------------------------------------
   Subroutine ReadPavr1x
      implicit none
      real(8):: alpha, r_e, fct
      integer:: i, j, k

      open(1,file="Data/Avr1x.dat",status='old')
         read(1,'(/(13E12.5))') (egg(k)        &
           , ((P(i,j,k),j=1,4),i=1,3),k=1,ne)
      close(1)
      print*, ' File Avr1x.dat is read'

      do k = 1,2
         lnegg(k) = dlog(egg(ne-2+k))
         do i = 1,nm
            do j = 1,nc
               lnP(i,j,k) = dlog(P(i,j,ne-2+k))
            enddo  ! j
         enddo  ! i
      enddo ! k

! Sigma_1_X.dat gives <(1-x)\sigma_{eep}>/(\alpha*r_e^2)
      open(1,file="Data/Sigma_1_X.dat",status='old')
         read(1,'(2E12.5)') (eep(k),avr1xeep(k),k=1,nee)
      close(1)
      print*, ' File Sigma_1_X.dat is read'
      alpha = 7.2973525698d-3	 ! $\alpha_{QED}$
      r_e   = 2.81794d-13        ! radius of electron in cm
      fct   = alpha*r_e**2*1d30  ! in $\mu$bn
      do k=1, nee
         eep     (k) = 1d-3*eep(k)     ! convert energy to GeV
         avr1xeep(k) = fct*avr1xeep(k) ! <1-x>\sigma(E) in $\mu$bn  
!		 write(*,'(2E12.5)') eep(k), avr1xeep(k)
      enddo
      lneep   (1) = dlog(eep     (nee-1))
      lneep   (2) = dlog(eep     (nee  ))
      lnavreep(1) = dlog(avr1xeep(nee-1))
      lnavreep(2) = dlog(avr1xeep(nee  ))

   end Subroutine ReadPavr1x

!-----------------------------------------

!-------------------------------------------------------------------------------
! Energy in GeV, <1-x>*\sigma in $\mu$bn for e^+ e^- photoproduction
   Subroutine Avr1xpair(Eg,avr)
      use msimsl, only: dcsiez
      implicit none
      real(8), intent( in):: Eg	
      real(8), intent(out):: avr	
      real(8),dimension(1):: xx, yy 
      real(8):: lneg  

      if(Eg.le.eep(1)) then
         avr = 0d0
      else if(Eg.le.eep(nee)) then
         xx(1) = Eg
         call dcsiez(nee,eep,avr1xeep,1,xx,yy)
         avr = yy(1)
      else
         lneg = dlog(Eg)
         avr  = dexp(FLIN(lneg,lneep(1),lneep(2)  &
              , lnavreep(1), lnavreep(2) ))
      endif

   end Subroutine Avr1xpair

!------------------------------------------------------------------------
! $ \int_xmin(E)^xmax(E) dx (1-x)^m d\sigma*E,x)/dx $  in $\m$bn
   real(8) Function Avr1xMI(Eg,im,ic)
      use gazlib, only: divdif
      implicit none
	  integer,intent(in):: im, ic  ! im = 1,2,3; ic =1,2 for pp and pn
      real(8),intent(in):: Eg		 ! energy in GeV
	  real(8),dimension(ne):: tmp

      if(Eg <= 0.15d0) then
         Avr1xMI = 0d0
      elseif(Eg > egg(ne)) then
         Avr1xMI  = dexp(FLIN(dlog(Eg),lnegg(1),lnegg(2)  &
                       , lnP(im,ic,1),lnP(im,ic,2)))
      else
	     tmp = P(im,ic,:) 
         Avr1xMI = DIVDIF(tmp,egg,ne,Eg,2)
      endif

   end Function Avr1xMI
!-------------------------------------------------------------------------------

!---------------------------------------------------
   real(8) Function FLIN(x,xl,xu,yl,yu)
      implicit none
      real(8), intent(in):: x, xl, xu, yl, yu

      FLIN = yl + (yu - yl)/(xu - xl)*(x - xl)

   end Function FLIN
!---------------------------------------------------

end module approx1x